##  Randomization inference analysis for EXPERIMENT 3 for sponsorship paper
rm(list = ls())

library(acs)
library(RCurl)
library(RJSONIO)


setwd("~/Dropbox/Current projects/kern-pietryka-crabtree/data/")
dat.a <- read.csv("sponsorship_3_March+15%2C+2018_08.01.csv", header = TRUE)
dat.a <- dat.a[-c(1:2),]

dim(dat.a)
colnames(dat.a)


##  drop test responses
dat.a <- dat.a[-c(1:15),]
dat.a$wave <- "A"


##  we need to add the second part of the survey
dat.b <- read.csv("sponsorship_3_v2_May+22%2C+2018_08.34.csv", header = TRUE)
dat.b <- dat.b[-c(1:2),]

dim(dat.b)
colnames(dat.b)


##  drop test responses
dat.b <- dat.b[-c(1:6),]
dat.b$wave <- "B"


##  which variable names have changed?
which(colnames(dat.a) != colnames(dat.b))


##  merge the two datasets, using colnames from the first survey
colnames(dat.b) <- colnames(dat.a)
dat <- rbind(dat.a,dat.b)
dim(dat)


##  consistency check: only real responses are included
table(dat$Status)


##  drop responses that were not completed
sel <- dat$Finished == "True"
table(sel)
dat <- dat[sel,]
##  In this study we have relatively large attrition, probably because of the
##  essay question



##  drop respondents who don't give consent
sel <- dat$Q1.2 == "Yes"
table(sel) 
dat <- dat[sel,]



##  one respondent answered twice for some reason. Drop second response
which(table(dat$mTurkCode) == 2)
sel <- which(dat$mTurkCode == "6600515")
dat <- dat[-sel[2],]

##  any respondents who answer multiple times?
stopifnot(length(dat$mTurkCode) == length(unique(dat$mTurkCode)))



##  drop respondents based on lat/long
latlong2fips <- function(latitude, longitude) {
  url <- "http://data.fcc.gov/api/block/find?format=json&latitude=%f&longitude=%f"
  url <- sprintf(url, latitude, longitude)
  json <- RCurl::getURL(url)
  json <- RJSONIO::fromJSON(json)
  as.character(json$County['FIPS'])
}
key <- api.key.install("f66712573be9fdfd1264a2d99bc0356a4b3edaf7")


###  this needs to be fixed to work properly; new web address

latlong2fips <- 
function(latitude, longitude) {
  url <- "https://geo.fcc.gov/api/census/block/find?latitude=%f&longitude=%f&format=json"
  url <- sprintf(url, latitude, longitude)
  json <- RCurl::getURL(url)
  json <- RJSONIO::fromJSON(json)
  as.character(json$County['FIPS'])
}


dat$state <- NA

m <- 1
for(m in 1:nrow(dat))	{

	temp <- latlong2fips(latitude = as.numeric(as.character(dat$LocationLatitude[m])),
						 longitude = as.numeric(as.character(dat$LocationLongitude[m])))

if(is.null(temp))	{	cat(m,as.numeric(as.character(dat$LocationLatitude[m])),
						as.numeric(as.character(dat$LocationLongitude[m])),"\n"); next() }
if(is.na(temp)) 	{	cat(m,as.numeric(as.character(dat$LocationLatitude[m])),
						as.numeric(as.character(dat$LocationLongitude[m])),"\n"); next() }
if(temp == "NULL")	{	cat(m,as.numeric(as.character(dat$LocationLatitude[m])),
						as.numeric(as.character(dat$LocationLongitude[m])),"\n"); next() }

	sel <- which(fips.state[1:51,1] == as.numeric(substr(temp,1,2)))
	temp <- fips.state[sel,2]

	if(length(temp) == 0)	{	cat(m,as.numeric(as.character(dat$LocationLatitude[m])),
								as.numeric(as.character(dat$LocationLongitude[m])),"\n");
								next() }
	dat$state[m] <- temp

	if(m %% 100 == 0) cat(m, "\n")

						}

##  drop respondents not from the 50 states
sel <- is.na(dat$state)
table(sel)
dat <- dat[!sel,]
dim(dat)

table(dat$state)


##  Deal with identical lat/long coordinates
##  Drop respondent if previous respondent with identical lat/long exists in data
dat$latlong <- paste(dat$LocationLatitude, dat$LocationLongitude, sep = "/")

dat$end.date <- strptime(dat[,2], format = "%Y-%m-%d %H:%M:%S")
dat <- dat[order(dat$end.date),]

#  consistency check
stopifnot(all(dat$end.date == sort(dat$end.date)))

dat$repeated <- NA
dat$repeated[1] <- FALSE


m <- 2
for(m in 2:nrow(dat))	{

	if(!any(dat$latlong[m] == dat$latlong[1:(m-1)])) { dat$repeated[m] <- FALSE } else
	{ dat$repeated[m] <- TRUE }

						}

table(dat$repeated)
sum(is.na(dat$repeated))
dat <- dat[!dat$repeated,]


##  drop responses recorded after Tillerson was fired on 13 March 2018 (8:44 AM - Mar 13, 2018)
dat <- dat[-c(1273:1366),]	# choose manually; make sure it stays consistent 



dat$know1 <- NA
dat$know2 <- NA
dat$know3 <- NA
dat$know4 <- NA
dat$know5 <- NA
dat$know6 <- NA

dat$know1.time <- NA
dat$know2.time <- NA
dat$know3.time <- NA
dat$know4.time <- NA
dat$know5.time <- NA
dat$know6.time <- NA

dat$Tr <- NA


m <- 1
for(m in 1:nrow(dat))	{

	# Senator term length
	temp <- 			c(	as.character(dat$Q3.1[m]),
							as.character(dat$Q11.2[m]),
							as.character(dat$Q17.2[m]))

	dat$know1[m] <- temp[temp != ""]
	dat$Tr[m] <- which(temp != "")	##  order: No logo, Harvard, FSU

	temp <- 			c(	as.character(dat$Q3.2_Page.Submit[m]),
							as.character(dat$Q11.3_Page.Submit[m]),
							as.character(dat$Q17.3_Page.Submit[m]))

	dat$know1.time[m] <- as.numeric(temp[temp != ""])


	# Secretary of State / Senator from Florida
	temp <- 			c(	as.character(dat$Q4.1[m]),
							as.character(dat$Q10.2[m]),
							as.character(dat$Q18.2[m]))

	dat$know2[m] <- temp[temp != ""]

	temp <- 			c(	as.character(dat$Q4.2_Page.Submit[m]),
							as.character(dat$Q10.3_Page.Submit[m]),
							as.character(dat$Q18.3_Page.Submit[m]))

	dat$know2.time[m] <- as.numeric(temp[temp != ""])


	# Federal spending
	temp <- 			c(	as.character(dat$Q5.1[m]),
							as.character(dat$Q12.2[m]),
							as.character(dat$Q19.2[m]))

	dat$know3[m] <- temp[temp != ""]

	temp <- 			c(	as.character(dat$Q5.2_Page.Submit[m]),
							as.character(dat$Q12.3_Page.Submit[m]),
							as.character(dat$Q19.3_Page.Submit[m]))

	dat$know3.time[m] <- as.numeric(temp[temp != ""])


	# Chief Justice
	temp <- 			c(	as.character(dat$Q6.1[m]),
							as.character(dat$Q13.2[m]),
							as.character(dat$Q20.2[m]))

	dat$know4[m] <- temp[temp != ""]

	temp <- 			c(	as.character(dat$Q6.2_Page.Submit[m]),
							as.character(dat$Q13.3_Page.Submit[m]),
							as.character(dat$Q20.3_Page.Submit[m]))

	dat$know4.time[m] <- as.numeric(temp[temp != ""])


	# Theresa May
	temp <- 			c(	as.character(dat$Q7.1[m]),
							as.character(dat$Q14.2[m]),
							as.character(dat$Q21.2[m]))

	dat$know5[m] <- temp[temp != ""]

	temp <- 			c(	as.character(dat$Q7.2_Page.Submit[m]),
							as.character(dat$Q14.3_Page.Submit[m]),
							as.character(dat$Q21.3_Page.Submit[m]))

	dat$know5.time[m] <- as.numeric(temp[temp != ""])


	# Attorney General
	temp <- 			c(	as.character(dat$Q8.1[m]),
							as.character(dat$Q15.2[m]),
							as.character(dat$Q22.2[m]))

	dat$know6[m] <- temp[temp != ""]

	temp <- 			c(	as.character(dat$Q8.2_Page.Submit[m]),
							as.character(dat$Q15.3_Page.Submit[m]),
							as.character(dat$Q22.3_Page.Submit[m]))

	dat$know6.time[m] <- as.numeric(temp[temp != ""])



	# Essay
	temp <- 			c(	as.character(dat$Q9.1[m]),
							as.character(dat$Q16.2[m]),
							as.character(dat$Q23.2[m]))

	dat$essay[m] <- temp[temp != ""]

	temp <- 			c(	as.character(dat$Q9.2_Page.Submit[m]),
							as.character(dat$Q16.3_Page.Submit[m]),
							as.character(dat$Q23.3_Page.Submit[m]))

	dat$essay.time[m] <- as.numeric(temp[temp != ""])

						}


table(dat$Tr)
##  test whether random assignment is as good as random
chisq.test(table(dat$Tr), simulate.p.value = TRUE, B = 10000)


table(dat$know1)
table(dat$know2)
table(dat$know3)
table(dat$know4)
table(dat$know5)
table(dat$know6)

hist(dat$know1.time)
hist(dat$know2.time)
hist(dat$know3.time)
hist(dat$know4.time)
hist(dat$know5.time)
hist(dat$know6.time)

hist(dat$essay.time)


##  Here we need to add distributions per treatment group as tables and plots
library(Hmisc)

out <- table(dat$know1, dat$Tr)
colnames(out) <- c("No logo","Harvard","Fitchburg State")
latex(out, file = "")


out <- table(dat$know2[dat$wave == "A"], dat$Tr[dat$wave == "A"])
colnames(out) <- c("No logo","Harvard","Fitchburg State")
latex(out, file = "")

out <- table(dat$know2[dat$wave == "B"], dat$Tr[dat$wave == "B"])
colnames(out) <- c("No logo","Harvard","Fitchburg State")
latex(out, file = "")


out <- table(dat$know3, dat$Tr)
colnames(out) <- c("No logo","Harvard","Fitchburg State")
latex(out, file = "")


out <- table(dat$know4, dat$Tr)
colnames(out) <- c("No logo","Harvard","Fitchburg State")
latex(out, file = "")


out <- table(dat$know5, dat$Tr)
colnames(out) <- c("No logo","Harvard","Fitchburg State")
latex(out, file = "")


out <- table(dat$know6, dat$Tr)
colnames(out) <- c("No logo","Harvard","Fitchburg State")
latex(out, file = "")



##  generate additional outcomes
trim <- function (x) gsub("^\\s+|\\s+$", "", x)

dat$know1 <- trim(dat$know1)
dat$know1.correct <- dat$know1 == "6 years"
table(dat$know1.correct)

dat$know1.dontknow <- dat$know1 == "Don't know"
table(dat$know1.dontknow)


dat$know2 <- trim(dat$know2)
dat$know2.correct <- dat$know2 == "Rex W. Tillerson" | dat$know2 == "Marco Rubio"
table(dat$know2.correct)

dat$know2.dontknow <- dat$know2 == "Don't know"
table(dat$know2.dontknow)


dat$know3 <- trim(dat$know3)
dat$know3.correct <- dat$know3 == "Foreign aid"
table(dat$know3.correct)

dat$know3.dontknow <- dat$know3 == "Don't know"
table(dat$know3.dontknow)


dat$know4 <- trim(dat$know4)
dat$know4.correct <- dat$know4 == "John Roberts"
table(dat$know4.correct)

dat$know4.dontknow <- dat$know4 == "Don’t know"
table(dat$know4.dontknow)


dat$know5 <- trim(dat$know5)
dat$know5.correct <- dat$know5 == "Prime Minister of the United Kingdom"
table(dat$know5.correct)

dat$know5.dontknow <- dat$know5 == "Don’t know"
table(dat$know5.dontknow)


dat$know6 <- trim(dat$know6)
dat$know6.correct <- dat$know6 == "Jeff Sessions"
table(dat$know6.correct)

dat$know6.dontknow <- dat$know6 == "Don't know"
table(dat$know6.dontknow)


##  consistency checks
stopifnot(sum(is.na(dat$know1.time)) == 0)
stopifnot(sum(is.na(dat$know2.time)) == 0)
stopifnot(sum(is.na(dat$know3.time)) == 0)
stopifnot(sum(is.na(dat$know4.time)) == 0)
stopifnot(sum(is.na(dat$know5.time)) == 0)
stopifnot(sum(is.na(dat$know6.time)) == 0)

stopifnot(sum(is.na(dat$essay.time)) == 0)


##  attentiveness
dat$attend <- trim(dat$Q24.1)
table(dat$attend)
out <- table(dat$attend, dat$Tr)
colnames(out) <- c("No logo","Harvard","Fitchburg State")
out <- rbind(out,NA)
rownames(out)[11] <- "% correct recall"
out[11,1] <- round(out[6,1] / sum(out[1:10,1]) * 100)
out[11,2] <- round(out[5,2] / sum(out[1:10,2]) * 100)
out[11,3] <- round(out[3,3] / sum(out[1:10,3]) * 100)

latex(out, file = "")



##  aggregate knowledge items
dat$dontknow <- dat$know1.dontknow +
				dat$know2.dontknow +
				dat$know3.dontknow +
				dat$know4.dontknow +
				dat$know5.dontknow +
				dat$know6.dontknow

dat$knowtime <- dat$know1.time +
				dat$know2.time +
				dat$know3.time +
				dat$know4.time +
				dat$know5.time +
				dat$know6.time

dat$knowtime.r <- rank(dat$knowtime)

dat$knowcorrect <-	dat$know1.correct +
					dat$know2.correct +
					dat$know3.correct +
					dat$know4.correct +
					dat$know5.correct +
					dat$know6.correct

table(dat$dontknow)
table(dat$knowcorrect)



##  outcomes by treatment group
out <- table(dat$dontknow, dat$Tr)
colnames(out) <- c("No logo","Harvard","Fitchburg State")
latex(out, file = "")

out <- table(dat$knowcorrect, dat$Tr)
colnames(out) <- c("No logo","Harvard","Fitchburg State")
latex(out, file = "")



##  Flesch scores and number of characters in essays
library(hunspell)
library(sylcount)

dat$char.no <- NA
dat$Flesch <- NA

m <- 1
for(m in 1:nrow(dat))	{

	##  Here we're adding a full stop because readability depends on full sentences.
	##  This does not affect scores if there already is a full stop.
	##  In contrast to pre-registration we do not check for proper English words.
	readability.out <- readability(paste(dat$essay[m],".", sep = ""))
	
	dat$char.no[m] <- unlist(readability.out[2])
	dat$Flesch[m] <- unlist(readability.out[8])

						}

stopifnot(sum(is.na(dat$char.no)) == 0)
hist(dat$char.no)

stopifnot(sum(is.na(dat$Flesch)) == 0)
range(dat$Flesch)
dat$Flesch[dat$Flesch == -Inf] <- max(dat$Flesch[dat$Flesch > -Inf])
hist(dat$Flesch)

# save.image("Experiment 3 intermediate results.RData")



##  randomization inference
##  compute sample test statistic for all six outcomes and all three-way comparisons
s <- 1000000


##  knowledge items
t.stat.dontknow.1v2 <- abs(mean(dat$dontknow[dat$Tr == 1]) - mean(dat$dontknow[dat$Tr == 2]))
t.stat.dontknow.1v3 <- abs(mean(dat$dontknow[dat$Tr == 1]) - mean(dat$dontknow[dat$Tr == 3]))
t.stat.dontknow.2v3 <- abs(mean(dat$dontknow[dat$Tr == 2]) - mean(dat$dontknow[dat$Tr == 3]))

t.stat.knowtime.1v2 <- abs(mean(dat$knowtime.r[dat$Tr == 1]) - mean(dat$knowtime.r[dat$Tr == 2]))
t.stat.knowtime.1v3 <- abs(mean(dat$knowtime.r[dat$Tr == 1]) - mean(dat$knowtime.r[dat$Tr == 3]))
t.stat.knowtime.2v3 <- abs(mean(dat$knowtime.r[dat$Tr == 2]) - mean(dat$knowtime.r[dat$Tr == 3]))

t.stat.knowcorrect.1v2 <- abs(mean(dat$knowcorrect[dat$Tr == 1]) - mean(dat$knowcorrect[dat$Tr == 2]))
t.stat.knowcorrect.1v3 <- abs(mean(dat$knowcorrect[dat$Tr == 1]) - mean(dat$knowcorrect[dat$Tr == 3]))
t.stat.knowcorrect.2v3 <- abs(mean(dat$knowcorrect[dat$Tr == 2]) - mean(dat$knowcorrect[dat$Tr == 3]))


##  essay question
dat$char.no.r <- rank(dat$char.no)
dat$essay.time.r <- rank(dat$essay.time)
dat$Flesch.r <- rank(dat$Flesch)

t.stat.char.no.1v2 <- abs(mean(dat$char.no.r[dat$Tr == 1]) - mean(dat$char.no.r[dat$Tr == 2]))
t.stat.char.no.1v3 <- abs(mean(dat$char.no.r[dat$Tr == 1]) - mean(dat$char.no.r[dat$Tr == 3]))
t.stat.char.no.2v3 <- abs(mean(dat$char.no.r[dat$Tr == 2]) - mean(dat$char.no.r[dat$Tr == 3]))

t.stat.essay.time.1v2 <- abs(mean(dat$essay.time.r[dat$Tr == 1]) - mean(dat$essay.time.r[dat$Tr == 2]))
t.stat.essay.time.1v3 <- abs(mean(dat$essay.time.r[dat$Tr == 1]) - mean(dat$essay.time.r[dat$Tr == 3]))
t.stat.essay.time.2v3 <- abs(mean(dat$essay.time.r[dat$Tr == 2]) - mean(dat$essay.time.r[dat$Tr == 3]))

t.stat.Flesch.1v2 <- abs(mean(dat$Flesch.r[dat$Tr == 1]) - mean(dat$Flesch.r[dat$Tr == 2]))
t.stat.Flesch.1v3 <- abs(mean(dat$Flesch.r[dat$Tr == 1]) - mean(dat$Flesch.r[dat$Tr == 3]))
t.stat.Flesch.2v3 <- abs(mean(dat$Flesch.r[dat$Tr == 2]) - mean(dat$Flesch.r[dat$Tr == 3]))


t.stat.rand <- matrix(NA,s,18)


set.seed(1)
m <- 1
for(m in 1:s)	{

	dat$Tr.r <- sample(c(1,2,3), nrow(dat), replace = TRUE)

	t.stat.rand[m,1] <- abs(mean(dat$dontknow[dat$Tr.r == 1]) - mean(dat$dontknow[dat$Tr.r == 2]))
	t.stat.rand[m,2] <- abs(mean(dat$dontknow[dat$Tr.r == 1]) - mean(dat$dontknow[dat$Tr.r == 3]))
	t.stat.rand[m,3] <- abs(mean(dat$dontknow[dat$Tr.r == 2]) - mean(dat$dontknow[dat$Tr.r == 3]))

	t.stat.rand[m,4] <- abs(mean(dat$knowtime.r[dat$Tr.r == 1]) - mean(dat$knowtime.r[dat$Tr.r == 2]))
	t.stat.rand[m,5] <- abs(mean(dat$knowtime.r[dat$Tr.r == 1]) - mean(dat$knowtime.r[dat$Tr.r == 3]))
	t.stat.rand[m,6] <- abs(mean(dat$knowtime.r[dat$Tr.r == 2]) - mean(dat$knowtime.r[dat$Tr.r == 3]))

	t.stat.rand[m,7] <- abs(mean(dat$knowcorrect[dat$Tr.r == 1]) - mean(dat$knowcorrect[dat$Tr.r == 2]))
	t.stat.rand[m,8] <- abs(mean(dat$knowcorrect[dat$Tr.r == 1]) - mean(dat$knowcorrect[dat$Tr.r == 3]))
	t.stat.rand[m,9] <- abs(mean(dat$knowcorrect[dat$Tr.r == 2]) - mean(dat$knowcorrect[dat$Tr.r == 3]))


	t.stat.rand[m,10] <- abs(mean(dat$char.no.r[dat$Tr.r == 1]) - mean(dat$char.no.r[dat$Tr.r == 2]))
	t.stat.rand[m,11] <- abs(mean(dat$char.no.r[dat$Tr.r == 1]) - mean(dat$char.no.r[dat$Tr.r == 3]))
	t.stat.rand[m,12] <- abs(mean(dat$char.no.r[dat$Tr.r == 2]) - mean(dat$char.no.r[dat$Tr.r == 3]))

	t.stat.rand[m,13] <- abs(mean(dat$essay.time.r[dat$Tr.r == 1]) - mean(dat$essay.time.r[dat$Tr.r == 2]))
	t.stat.rand[m,14] <- abs(mean(dat$essay.time.r[dat$Tr.r == 1]) - mean(dat$essay.time.r[dat$Tr.r == 3]))
	t.stat.rand[m,15] <- abs(mean(dat$essay.time.r[dat$Tr.r == 2]) - mean(dat$essay.time.r[dat$Tr.r == 3]))

	t.stat.rand[m,16] <- abs(mean(dat$Flesch.r[dat$Tr.r == 1]) - mean(dat$Flesch.r[dat$Tr.r == 2]))
	t.stat.rand[m,17] <- abs(mean(dat$Flesch.r[dat$Tr.r == 1]) - mean(dat$Flesch.r[dat$Tr.r == 3]))
	t.stat.rand[m,18] <- abs(mean(dat$Flesch.r[dat$Tr.r == 2]) - mean(dat$Flesch.r[dat$Tr.r == 3]))


	if(m %% 1000 == 0)	{ cat(m, "\n") }

				}

# save.image("Experiment 3 results.RData")
# load("Experiment 3 results.RData")



p <- matrix(NA,6,4)
rownames(p) <- c("% Don't know","Response time","% Correct answers","Number of characters","Response time","Flesch score")
colnames(p) <- c("No logo vs. Harvard","No logo vs. FSU","FSU vs. Harvard","Joint")

##  p-values for each outcome and each two-way comparison
p[1,1:3] <- 

c(
mean(t.stat.dontknow.1v2 < t.stat.rand[,1]),
mean(t.stat.dontknow.1v3 < t.stat.rand[,2]),
mean(t.stat.dontknow.2v3 < t.stat.rand[,3]))

p[2,1:3] <- 

c(
mean(t.stat.knowtime.1v2 < t.stat.rand[,4]),
mean(t.stat.knowtime.1v3 < t.stat.rand[,5]),
mean(t.stat.knowtime.2v3 < t.stat.rand[,6]))

p[3,1:3] <- 

c(
mean(t.stat.knowcorrect.1v2 < t.stat.rand[,7]),
mean(t.stat.knowcorrect.1v3 < t.stat.rand[,8]),
mean(t.stat.knowcorrect.2v3 < t.stat.rand[,9]))



p[4,1:3] <- 

c(
mean(t.stat.char.no.1v2 < t.stat.rand[,10]),
mean(t.stat.char.no.1v3 < t.stat.rand[,11]),
mean(t.stat.char.no.2v3 < t.stat.rand[,12]))

p[5,1:3] <- 

c(
mean(t.stat.essay.time.1v2 < t.stat.rand[,13]),
mean(t.stat.essay.time.1v3 < t.stat.rand[,14]),
mean(t.stat.essay.time.2v3 < t.stat.rand[,15]))

p[6,1:3] <- 

c(
mean(t.stat.Flesch.1v2 < t.stat.rand[,16]),
mean(t.stat.Flesch.1v3 < t.stat.rand[,17]),
mean(t.stat.Flesch.2v3 < t.stat.rand[,18]))


##  Using Young's approach to combine test statistics
Omega.1 <- solve(var(t.stat.rand[,1:3]))
Omega.2 <- solve(var(t.stat.rand[,4:6]))
Omega.3 <- solve(var(t.stat.rand[,7:9]))
Omega.4 <- solve(var(t.stat.rand[,10:12]))
Omega.5 <- solve(var(t.stat.rand[,13:15]))
Omega.6 <- solve(var(t.stat.rand[,16:18]))


t.stat.1 <-	t(c(t.stat.dontknow.1v2,t.stat.dontknow.1v3,t.stat.dontknow.2v3)^2) %*%
			Omega.1 %*%
		  t(t(c(t.stat.dontknow.1v2,t.stat.dontknow.1v3,t.stat.dontknow.2v3)^2))

t.stat.2 <-	t(c(t.stat.knowtime.1v2,t.stat.knowtime.1v3,t.stat.knowtime.2v3)^2) %*%
			Omega.2 %*%
		  t(t(c(t.stat.knowtime.1v2,t.stat.knowtime.1v3,t.stat.knowtime.2v3)^2))

t.stat.3 <-	t(c(t.stat.knowcorrect.1v2,t.stat.knowcorrect.1v3,t.stat.knowcorrect.2v3)^2) %*%
			Omega.3 %*%
		  t(t(c(t.stat.knowcorrect.1v2,t.stat.knowcorrect.1v3,t.stat.knowcorrect.2v3)^2))

t.stat.4 <-	t(c(t.stat.char.no.1v2,t.stat.char.no.1v3,t.stat.char.no.2v3)^2) %*%
			Omega.4 %*%
		  t(t(c(t.stat.char.no.1v2,t.stat.char.no.1v3,t.stat.char.no.2v3)^2))

t.stat.5 <-	t(c(t.stat.essay.time.1v2,t.stat.essay.time.1v3,t.stat.essay.time.2v3)^2) %*%
			Omega.5 %*%
		  t(t(c(t.stat.essay.time.1v2,t.stat.essay.time.1v3,t.stat.essay.time.2v3)^2))

t.stat.6 <-	t(c(t.stat.Flesch.1v2,t.stat.Flesch.1v3,t.stat.Flesch.2v3)^2) %*%
			Omega.6 %*%
		  t(t(c(t.stat.Flesch.1v2,t.stat.Flesch.1v3,t.stat.Flesch.2v3)^2))


out.1 <- rep(NA,nrow(t.stat.rand))
out.2 <- rep(NA,nrow(t.stat.rand))
out.3 <- rep(NA,nrow(t.stat.rand))
out.4 <- rep(NA,nrow(t.stat.rand))
out.5 <- rep(NA,nrow(t.stat.rand))
out.6 <- rep(NA,nrow(t.stat.rand))

m <- 1
for(m in 1:nrow(t.stat.rand))	{

	out.1[m] <- t.stat.rand[m,1:3]^2   %*% Omega.1 %*% t(t(t.stat.rand[m,1:3]^2))
	out.2[m] <- t.stat.rand[m,4:6]^2   %*% Omega.2 %*% t(t(t.stat.rand[m,4:6]^2))
	out.3[m] <- t.stat.rand[m,7:9]^2   %*% Omega.3 %*% t(t(t.stat.rand[m,7:9]^2))
	out.4[m] <- t.stat.rand[m,10:12]^2 %*% Omega.4 %*% t(t(t.stat.rand[m,10:12]^2))
	out.5[m] <- t.stat.rand[m,13:15]^2 %*% Omega.5 %*% t(t(t.stat.rand[m,13:15]^2))
	out.6[m] <- t.stat.rand[m,16:18]^2 %*% Omega.6 %*% t(t(t.stat.rand[m,16:18]^2))

								}

p[1,4] <- mean(as.vector(t.stat.1) < out.1)
p[2,4] <- mean(as.vector(t.stat.2) < out.2)
p[3,4] <- mean(as.vector(t.stat.3) < out.3)
p[4,4] <- mean(as.vector(t.stat.4) < out.4)
p[5,4] <- mean(as.vector(t.stat.5) < out.5)
p[6,4] <- mean(as.vector(t.stat.6) < out.6)


p2 <- round(p,2)

library(Hmisc)
latex(p2, file = "")



##  Omnibus test using Young's approach
Omega <- solve(var(t.stat.rand))

t.stat <-	t(c(t.stat.dontknow.1v2,t.stat.dontknow.1v3,t.stat.dontknow.2v3,
				t.stat.knowtime.1v2,t.stat.knowtime.1v3,t.stat.knowtime.2v3,
				t.stat.knowcorrect.1v2,t.stat.knowcorrect.1v3,t.stat.knowcorrect.2v3,
				t.stat.char.no.1v2,t.stat.char.no.1v3,t.stat.char.no.2v3,
				t.stat.essay.time.1v2,t.stat.essay.time.1v3,t.stat.essay.time.2v3,
				t.stat.Flesch.1v2,t.stat.Flesch.1v3,t.stat.Flesch.2v3)^2) %*%
				Omega %*%
		  t(t(c(t.stat.dontknow.1v2,t.stat.dontknow.1v3,t.stat.dontknow.2v3,
				t.stat.knowtime.1v2,t.stat.knowtime.1v3,t.stat.knowtime.2v3,
				t.stat.knowcorrect.1v2,t.stat.knowcorrect.1v3,t.stat.knowcorrect.2v3,
				t.stat.char.no.1v2,t.stat.char.no.1v3,t.stat.char.no.2v3,
				t.stat.essay.time.1v2,t.stat.essay.time.1v3,t.stat.essay.time.2v3,
				t.stat.Flesch.1v2,t.stat.Flesch.1v3,t.stat.Flesch.2v3)^2))

out <- rep(NA,nrow(t.stat.rand))

m <- 1
for(m in 1:nrow(t.stat.rand))	{

	out[m] <- t.stat.rand[m,]^2 %*% Omega %*% t(t(t.stat.rand[m,]^2))

								}

round(mean(as.vector(t.stat) < out),2)











.
